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The study of metal nanoparticles plays a central role in the emerging novel technologies employing optics 
beyond the diffraction limit. Combining strong surface plasmon resonances, high intrinsic nonlinearities 
and deeply subwavelength scales, arrays of metal nanoparticles offer a unique playground to develop novel 
concepts for light manipulation at the nanoscale. Here we suggest a novel principle to control localized 
optical energy in chains of nonlinear subwavelength metal nanoparticles based on the fundamental 
nonlinear phenomenon of modulation instability. In particular, we demonstrate that modulation instability 
can lead to the formation of long-lived standing and moving nonlinear localized modes of several distinct 
types such as bright and dark solitons, oscillons, and domain walls. We analyze the properties of these 
nonlinear localized modes and reveal different scenarios of their dynamics including transformation of one 
type of mode to another. We believe this work paves a way towards the development of nonlinear 
nanophotonics circuitry. 



Modulation instability (MI) is a nonlinearity-induced phenomenon that can be observed in many different 
branches of physics^ Due to MI, small noise-driven amplitude and phase perturbations rapidly boost 
under the influence of nonlinearity and diffraction (or dispersion). Consequently, broad optical beams 
(or quasicontinuous wave pulses) decay into optical filaments (or pulse trains)^ MI also shows up in the pattern 
formation from noise for partially spatially incoherent light beams propagating through noninstantaneous non- 
linear media^. Finally, MI development in granular and liquid media is often associated with generation of 
nonlinear localized states in the form of oscillons (or oscillating solitons) and domain walls^"^. Hopefully, 
consideration of such intriguing phenomena in plasmonic systems such as metallic nanowires or nanoparticle 
arrays may open novel rotes for light control at the deeply subwavelength structures. 

Owing to the ability to squeeze light into a nanometer size regime, surface plasmons supported by nobel 
metal nanostructures are exploited in a huge number of applications ranging from photonic nanocircuitry to 
biological sensors However, another remarkable property of nobel metals - extremely high intrinsic 
nonlinearity - remained obscure for a long time because of their natural optical opaqueness. Following in 
the footsteps of ground-breaking works on a cubic susceptibility of metal colloids^^"^^, recent advances in the 
development of high-quality nanostructured systems manifest how to engineer the field penetration inside the 
nanostructure to fully use its nonlinear response^^. In this context, one may mention optical limiting and self- 
phase modulation in arrays of structured nanoparticles^^, second and third harmonic generation in nanos- 
tructured metal films and nanoantennas^^"^^ subwavelength solitons in metaldielectric multilayers^^"^^ and 
arrays of metal nanowires^^ plasmonic kinks and oscillons in nanoparticle arrays^^'^^, ultrafast all-optical 
switching in metallic photonic crystals^^ as well as metal-semiconductor and metal -dielectric nanoanten- 
nas^^~^\ generation of THz radiation from silver nanodimers^^'^^, and a nanoradar with periodically rotating 
scattering pattern^^. 

Here we present a comprehensive study of MI in subwavelength nonlinear systems for a chain of optically 
driven metallic nanoparticles with a Kerr-like nonlinear response. We demonstrate the existence of novel types of 
nonlinear effects in such systems, including generation of domain walls (or switching waves) as well as long-lived 
standing and moving nonlinear localized modes in the form of bright and dark oscillons and solitons. We reveal 
that a wide variety of scenarios of MI development allows for the mode transformation from one type to another. 
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changing the propagation direction and a velocity of drifting oscil- 
lons, sohtons, and domain walls as well as the formation of stable 
domain walls connecting not only different stationary states but also 
the states with different types of nonlinear dynamics. 

Results 

Model and governing equations. We consider a chain of identical 
spherical nanoparticles placed close to each other and embedded into 
a Si02 host medium with a permittivity Sh = 2.15 (see Fig. 1). We 
assume that the particle radius and the center-to- center spacing are 
a = 10 nm and d = 30 nm, respectively. Ratio a/d satisfies the 
condition a/d < 1/3, so that we can employ the point dipole 
approximation^^. Assuming the nanoparticles made of silver with a 
nonlinear Kerr-like response, we take the dielectric constant in the 

form ^Ag ^^Ag + X*^^^ I , where the linear part is given by the 
Drude formula, s\^ = £oo [ft^(cL> — /v)], with goo = 4.96, ^cl)^= 

9.54 eV, and 0.055 eV'^ and E^^^^^ is the local field inside the nth 
particle. In general, the value of the cubic susceptibility for metallic 
nanoparticles dependents on a number of factors, including the type 
of metal, particle size, external pulse duration and frequency, and 
others^^. However, the analytical quantum model derived in Refs.^^'^^ 
and confirmed with numerical simulations^^ and experimental 
j^^^i2,i4 si^owed that silver nanoparticles with radius 10 nm and 
driven at the frequency close to the frequency of the surface 
plasmon resonance possess a remarkably high and purely real 

cubic susceptibility j^^^ ^3 x 10 ~^ esu. Note fused silica possesses 
negligibly weak cubic nonlinearity (~ 10"^^ esu^^) as well as a good 
enough optical transparency^^ to maintain the strong laser powers 
needed for the observation of MI. 

We study a nanoparticle chain driven by an arbitrary optical field 
with the frequency close to the frequency of the surface plasmon 



I 




Particle number 

Figure 1 | Schematic of a chain of metallic nanoparticles illuminated by a laser beam and profiles of nonlinear localized states, (a) Arrows indicate 
particle polarizations for a bright soliton/oscillon. Panels (b), (c) and (d) depict profiles of the polarizations for a typical soliton/oscillon configuration of 
bright and dark forms as well as domain wall, respectively. 



resonance of an individual particle, and analyze the dynamical response 
of the particle polarizations, p„, through the model proposed in our 
previous article^^, which yields the following system of coupled equa- 
tions for the slowly varying amplitudes of the particle dipole moments. 



Ar>l- 



dT 

where 



/ ^ + ( - /y + a + |P„ |^)P« + ^ Gl^Pl = 4, 



(1) 



. _ / ikod 1 A e-*"^!^-'"! 

P^'" and are dimensionless slowly varying amplitudes of the 
particle dipole moments and external electric field, respectively, the 
indices '±' and 'H' stand for the transverse and longitudinal compo- 
nents with respect to the chain axis, r| defines strength of dipole-dipole 

interaction, |P„|^ = |P;J- 1^ + |p|| |^ is a nonlinear term appeared after 
expressing E^^^^^ via pn, y describes thermal and radiation losses of 
particles, ko = coo/ c^/sh is the wavenumber of light in a host matrix, 
Q = (co - (Do)/oDo, and x = a^ot (see Methods). We notice that this 
model involves all particle interactions through the dipole fields, and it 
can be applied to both fmite and infmite chains. 

To begin with, we excite the infmite chain by an electric field with 
one of the two polarizations: (i) E„ = [Eq exp{ — iQdn),0) and (ii) 
En= [0,Eo exp( — /Q(in) J, where Q is the wavenumber along the 
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chain axis. We suppose that Q may vary in a wide range, 0 < Q(i < 7i. 
In practice 0 < Q < /cq can be reaUzed through tilted Ught incidence 
on the chain, while Q > ko can be provided by means of a scheme 
based on, e.g., total internal reflection, as it was suggested recently for 
investigation of localized plasmon modes in metallic nanoparticle 
chains^^. The simplest case of homogenous excitation, i.e. Q = 0, 
has been considered earlier^^, and we observed periodic pattern 
formation and generation of oscillons. Below we will show that 
adjustment of Q opens a wide variety of novel scenarios of MI 
development. 

Spatially modulated external field allows us to treat particle dipole 
moments as P;['" =-Pj"'" exp{ — iQdn). In this case, the system sta- 
tionary states can be written as follows 



-ly 



pi,ll_£i,l 



(2) 



where 



and 



(kodf 



ikod 



cos ( Qdf) exp ( — //co df) , 



ikod 1 



cos{Qdj)exp{ — ikodj) . 



A transition from G^jj^ to A^^'" has been made via the replacement | n 
— m\ = j and taking into account symmetry structure of the series. 
This equation has multiple solutions for Q< — Re ^j^^ A^^'" — 

|y — Im Y^J^ J A^^'" I leading to a bistable regime. 

Linear stability analysis and numerical simulations of modulation 
instability. Next, we analyze linear stability of the stationary states 
(2) with respect to weak spatiotemporal perturbations, and derive the 
expression for the instability growth rate (see Methods), 



1=^1,11 



2 Pit 



' + Q + Re^5/' 



1/2 



where y^ ||=Im J^^^ 



1^; 



^,11 



-y, B:^'"=A:^'" cos{Kdj). Thus, the 



initial nonlinear modulated states (2) become unstable provided 

7-1,11 



> 0. The stability depends on the external field parameters 
Q, and Q as well as on the modulation wavenumber of 



perturbation, K. 

Consider the case of the transversal excitation in detail. Figure 2(a) 
shows regions of MI in the plane (Q, \Eq\ ) obtained from the 
condition lj_ = 0 at any K and Qd= 0.7 that corresponds to grazing 
incidence of light at the angle 2.5° relative to the chain axis. 
Importantly, external field modulation results in the remarkable fact 
that the MI development can be reached with respect to both fast 
(with K < ko) and slow (with K > ko) system eigenmodes in contrast 
to the case of homogeneous excitation, for which MI generation is 
achievable with respect to slow eigenmodes only^^. The first possibil- 
ity opens a promising perspective for dynamical engineering of the 
chain scattering pattern similar to nanoradar^^; whereas the second 
provides the straightforward rote for manipulation of strongly loca- 
Uzed nonlinear states. In this work we focus on the latter approach. 

Contour maps of lj_ in the planes {Kd, \Eq\ ) and {Qd, Kd) are 
presented in Figs. 2(b) and 2(c), respectively. As follows from these 
figures, one can generate MI within either slow or fast eigenmodes by 
proper adjustment of | and Q, enabling dynamical transforma- 
tion of free propagating radiation into a nonlinear localized form and 
vice versa, a localized external excitation into free propagating radi- 
ation. Theses possibilities correspond to the zones with K> ko,Q< 
ko and K < ko, Q> ko in Fig. 2(c), respectively. 
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Figure 2 | Bifurcation diagram and contour maps of the instability 
growth rate at the transversal excitation, (a) Bifurcation diagram on the 
parameter plane of Q = (co — coq)/(jOo and \Eq | and Qd = 0.7, with the 
regions of bistability (blue) and modulation instability (green and red). 
Red and green zones correspond to MI development relative to slow (with 
K > ko) and fast (with K < ko) eigenmodes of the chain, respectively. The 
purple star denotes intensity \Eq\ = 10~^ (which is also marked by the 
horizontal dashed line in (b)) and frequency Q = —0.03 used for drawing 
contour maps of /l^ on the planes (b) {Kd,\EQ | ) and (c) {Kd, Qd) as well 
as in numerical simulations of Eq. (1) shown in Fig. 3. Inset in (b) 
demonstrates dependence of the polarization \Po\^ on |^ for Q = 
— 0.03. Red and green dots indicate red and green zones of modulation 
instability from (a), respectively. Inset in (c) shows the enlarged contour 
map of 2^ in the vicinity of Qd = 0.7. 

In order to examine the spatiotemporal evolution of the chain 
beyond the instability point, we perform numerical simulations of 
Eq. (1) for a finite chain (with 100 nanoparticles) at zero initial 
conditions, when the external field increases slowly approaching to 
the saturation level Eq, which is reached at t ~ 100. We selected 
= 10-5, Q = -0.03 and Q = 0.7 to undergo MI out of the 
bistability zone and with respect to slow eigenmodes only (see 
Figs. 2(a), (b), and (c)). Edge effects acted as small perturbations 
needed for MI onset. 

The result is presented in Fig. 3, where we observe that instability 
provokes formation of a train of bright oscillons. Oscillons appear at 
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Figure 3 | Train of bright oscillons. Dynamics of the particle polarizations 

obtained by numerical simulations of Eq. (1) at Q = —0.03, 
Eq = 10~^, and Qd = 0.7, when the light beam propagating at the angle 
2.5° with respect to the chain axis provokes generation of a train of bright 
oscillons. Panel (a) shows a snapshot of |P;|- 1 at t = 5000. Panels (b) and 
(c) represent 2D and 3D views of the spatiotemporal dynamics of | , 
respectively. See also Supplementary Material for the associated time 
animations of particle polarizations. 



the right side of the chain and move with the speed ~ 0.013c in the 
opposite direction with respect to the direction of light propagation, 
as shown in Fig. 1. Having reached the left side of the chain, they 
transform into free propagating radiation, and a small part of energy 
is reflected back in the chain, eventually decaying. We notice this type 
of oscillons is generated in a two-pick form, and the first oscillon 
turned into an oscillon and a soliton, which later transformed back 
into an oscillon. Remarkably, all bright oscillons are localized on 7 
particles, and their width is about 210 nm. 

Next, we conduct the similar analysis for the case of the longit- 
udinal excitation. Figures 4(a) to (c) show the corresponding bifurca- 
tion diagram and contour maps. In contrast to the transversal case, 
the MI growth of both localized (with Q > /cq) and free propagating 
(with Q < /cq) excitations can be reached within slow eigenmodes 
only. Since the instability growth rate takes the maximum values for 
modes with Kd ~ 7i, one may expect that the subsequent evolution of 
the unstable system will lead to formation of a high-order mode. 

Indeed, numerical simulations confirm appearance of small- 
amplitude high-order modulations, but accompanied by an array 
of standing dark sohtons, as shown in Figs. 5 (a) and (b). 
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Figure 4 | Bifurcation diagram and contour maps of the instability 
growth rate at the longitudinal excitation, (a) Bifurcation diagram at Qd 

I II P 

= 0.7. The purple star denotes intensity £q = 5 x 10 ~^ (which is also 

marked by the horizontal dashed line in (b)) and frequency Q = —0.06 

I II P 

used for drawing contour maps of 2 on the planes (b) {Kd, LEq ) 

(c) {Kd, Qd) as well as in numerical simulations of Eq. (1) shown in Figs. 5 

I II P 

(a) and (b). Inset demonstrates dependence of the polarization \PU on 
1^0 1 for Q = -0.06. Red dots indicate the area of MI. 

Surprisingly, despite long-range interactions, all these solitons are 
localized only on 1 particle. Nevertheless, nonlocality of particle 
interaction is manifested in a soliton-to-soliton distance which varies 
between 8 and 12 particles in the middle of the chain and decreases to 
2-6 nanoparticles close to the edges. 

In addition to dark solitons, we observed generation of standing 
dark oscillons, as shown in Figs. 6(a) and (b). They also are locahzed 
on 1 particle, but with a flat background. Although MI growth leads 
to formation of dark oscillons along all the chain, only 3 of them 
eventually remain stable; while others decay. Therefore, one may 
conclude that in general standing dark oscillons are less stable than 
standing dark solitons. 

Figures 6 (c) and (d) represent a train of dark oscillons. Similar to 
drifting bright oscillons, drifting dark oscillons are generated at the 
right edge of the chain and move with the starting velocity 0.014c, 
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Figure 5 | Dark solitons observed in the longitudinally excited chain at fl 
= -0.06, I = 5 X 10"^, Qd = 0.7. Panel (a) shows a snapshot of the 

particle polarizations \P\^ at t = 5000; while panel (b) represents a 2D 
view of the spatiotemporal dynamics of |pU | . See also Supplementary 
Material for the associated time animation. 



slowly decelerating. Then they collapse or stay (not shown) on the 
way, not reaching the other side. Interestingly, this kind of oscillons 
combines a narrow width about 3 particles (apart from the first 
oscillon whose width varies from 1 1 to 3 particles) with quite a long 
tale consisting 1 1 particles. 

The important property of oscillons is the ability to change the 
movement direction. That case is presented in Figs. 6 (e) and (f), 
where the walking oscillon in a high-contrast background is slow 
drifting along the chain, switching the movement direction from 
time to time. Remarkably, the absolute value of the velocity remains 
constant and equals to 0.0037c. While this simulation shows only the 
oscillon which walks around the middle of the chain, other scenarios 
of MI development with close parameters allow arrest of the walking 
oscillon by the edge and subsequent transformation into the surface 
locahzed state (not shown). 

Domain walls. Finally, we analyze domain walls (or switching 
waves). In general there are two mechanisms leading to domain 
wall formation: bistability and modulation instability. Bistability- 
induced domain walls also often referred as kinks connect two 
different stable states of particle polarizations. MI makes twofold 
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Figure 6 | Examples of nonlinear localized states. Panels (a,b) indicate generation of dark oscillons in the longitudinally excited chain at Q = 

I l|2 2 

Eq =4.5 X 10^^, Qd = 2. Panels (c,d) show a train of dark oscillons in the transversally excited chain at Q = —0.06, \Eq | = 1.2 x 10~^, Qd 

Panels (e,f) represent the case of the oscillon with a changeable direction of movement but the constant absolute value of the speed ~ 0.0037c in the 
transversally excited chain at Q = —0.0475, =7 x 10~^,Q(i= 7i, respectively. Panels (g,h) demonstrate formation ofthe domain wall between parts 
of the chain with regular and chaotic-like dynamics. The horizontal dashed line in (h) marks a time moment for the frozen frame shown in (g). See also 
Supplementary Material for the associated time animations. 
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impact with respect to domain walls. When MI occupies at least one 
of the bistability branches, domain walls generation typically cannot 
be realized. For instance, this case corresponds to the bifurcation 
diagrams shown in Figs. 2 (a) and 4 (a). On the other hand, MI is 
able to generate domain walls beyond the bistability regime on its 
own. Figures 6 (g) and (h) demonstrate an example of such a domain 
wall which appeared between the states with high- order 
modulations, accompanied by small amplitude beatings, and a 
chaotic-like behavior. Having appeared at the left edge, such a 
switching wave drifts with a variable velocity until stay close to the 
other side. Interestingly, domain wall stopping results in changing 
chaotic-like dynamics with regular beatings. 

In our previous work^^ we showed that in the case of homogenous 
transversal excitation MI occupies only a small region of intensities 
close to the upper threshold of the bistability area, making possible 
kink generation. For that case we also observed kinks with zero 
velocity and one-way movement^^. Here we extend our analysis fur- 
ther and reveal switching the kink velocity from negative to positive 
values. The velocity sign is defined positive if the motion of the 
switching wave leads to the expansion of the region occupied by 
the upper branch of bistability; otherwise, the velocity is negative. 

To observe kink dynamics, we perform numerical simulations of 
Eq. (1) for a finite array (with 200 nanoparticles) at initial conditions 
corresponding to different branches of the homogeneous stationary 
solution (2) over different parts of the array. Light intensity is sup- 
posed to be a step-like function of time, as shown in Fig. 7 (a). 
Constant Eq results in the constant kink velocity, v. That is why a 
step-like time dependency of Eq allows us to control v. 

The results are presented in Fig. 7 (b) where the kink velocity is 
switched from negative to positive values, passing through zero. A 
rigorous analysis manifests that there is no smooth transition 
between kinks with negative, positive, and zero velocities, as shown 
in Fig. 7 (c). Moreover, these states possess different structures: kinks 
with positive and zero velocities demonstrate zero width; while the 
width of kinks with negative velocity extends for 12 particles. We also 
notice that negative -velocity kinks exist for the same intensities as 



stationary (zero -velocity) kinks do. Therefore, to provoke domain 
wall moving in the negative direction one should force the dipole 
next to the domain wall into the lower branch of the bistability band 
until the domain wall starts moving. An opposite transition can be 
realized through only variation in the intensity, as shown in Fig. 7 (c). 

Discussion 

The above examples demonstrate how modulation instability and 
bistability can be used to manipulate light in arrays of nonlinear plas- 
monic nanoparticles coupled through all dipole fields. In particular, 
we revealed the generation of long-lived standing and moving non- 
linear localized modes in the form of plasmonic oscillons and soli- 
tons of both bright and dark types as well as plasmonic domain walls. 

Notably, formation of solitons/oscillons can be predominantly 
governed by either short-range or combination of short-range and 
long-range interactions. This is best illustrated in localization ran- 
ging from 1 to 11 nanoparticles that answers the dimensional scales 
from 0.05/1 to 0.8>^, because Q ^ 1 corresponds to the radiation 
wavelength 1 ~ 400 nm. 

Long-range interaction also features in kink dynamics. Addressing 
to a similar model applied to investigation of kinks in arrays of split- 
ring resonators coupled only through near field^^'^^, we find that 
long-range interaction splits branches with positive and zero velo- 
cities for V ^ and extends the range of intensities for stationary 
kinks to the lower threshold of the bistability area. In addition, a 
near-field model does not distinguish the structure of the kinks with 
positive and negative velocities; while our model taking into account 
coupling via all dipole fields indicates a dramatic change in kink 
width when the velocity sign is changed (see above). 

It may be of interest to compare MI -induced plasmon solitons/ 
oscillons in nanoparticle chains with plasmon solitons supported by 
metal-dielectric multilayered structures'^ and arrays of metallic 
nanowires'^. All these nonlinear localized states demonstrate subwa- 
velength localization. However, solitons in ID multilayers and 2D 
nanowire lattices extend for 031 and 0.61, respectively; whereas the 
width of solitons/oscillons in nanoparticle chains may reach 0.05L 




Bistability area 



(b) 






lE^PxIO-^ 

Figure 7 | Control over the kink velocity. Panel (a) demonstrates stimuli vs. time resulting in spatiotemporal dynamics ofa kink (switching wave), shown 
in (b) (Q = — 0.1, Q = 0). Panel (c) indicates the normalized kink velocity vs. the intensity of the applied field. A blue curve does not reach the upper 
threshold of the bistability area because kinks do not exist here owing to ML Insets show kink structures with positive and negative velocities. See also 
Supplementary Material for the associated time animation. 
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Besides, MI may create dark solitons/oscillons which were never 
studied before for plasmonic systems. Thus, MI in arrays of plasmo- 
nic nanoparticles opens a perspective way for generation of extre- 
mely squeezed nonlinear localized states unreachable with other 
nonlinear plasmonic structures. 

It is also insightful to estimate the maximal pulse duration, because 

the saturation intensity of the external field lE'^ 



10- 



10- 



corresponded to the intensity of ~ 1 - 10 MW/cm^, can lead to 
thermal damage. Using the value of the ablation threshold of 
3.96 J/cm^ obtained for silver particles in a Si02 host matrix in the 
picosecond regime of illumination^^ and taking into account the 
amplification of the electric field inside the Ag nanoparticle due to 
surface plasmon resonance, we evaluate the maximal pulse duration 
as 5 ns for the intensity 0.75 MW/cm^ (|-^o" I = required for 

generation of a train of bright oscillons. Since a typical interval 
between formation of oscillons is about 150 fs (At ^700), one may 
observe generation of ~ 33000 oscillons before a chain will be 
destroyed. Figures 5 to 7 show similar characteristic times for other 
nonlinear localized states and domain walls. Therefore, all predicted 
phenomena can be readily observed in experiment. 

Our findings provide a starting platform for further experimental 
studies of the nonlinearity- induced instabilities and associated phe- 
nomena in plasmonic nanostructures which could have important 
implications for active nanophotonic devices operating beyond the 
diffraction limit. 

Methods 

In our theoretical analysis, we use following dimensionless quantities 



Soo -\-2sh) 



38h 



+ (M'- 



To derive the instability growth rate, we utilized the standard technique"*^, 
and take the small perturbations in the form of chain eigenmodes: 
(5P^'II =Ai I exp( — + t) u exp (^iK dn ^ ||T^, where means 
complex conjugation, K is the modulation wavenumber, Aj_^\\ are amplitudes of 
perturbation. Having substituted P^'' =P(i"'' exp( — zQc/n) + ^P^'" into Eq. (1), we 
come to the expression for the instability growth rate. 

Numerical simulations of Eq. (1) were performed by means of a commercial 
package Matlab v. 7 with a fourth-order Runge-Kutta scheme, and relative and 
absolute error tolerances were fixed as 10"^ and 10"^, respectively. 
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